Oyster exposure experiments, length and weight analyses

knitr::opts_chunk$set(echo = TRUE, warning = F, message = F, fig.path = 'figs/', dev.args = list(family = 'serif'))

library(tidyverse)
library(janitor)
library(readxl)
library(patchwork)
library(survival)
library(ggfortify)

source("R/funcs.R")

data(allind)
data(indbymods)
data(dirindbymods)

# consistent treatment terminology names
trts <- tibble(
  shrtlab = c("7.7C", "7.7A0.2", "7.7A0.5", "8.0C", "8.0A0.2"),
  lngslab = c( "7.7 Constant", "7.7 Fluctuating 0.2A", "7.7 Fluctuating 0.5A", "8.0 Constant", "8.0 Fluctuating 0.2A")
)

# week combinations
wks <- sort(unique(allind$week)) %>% 
  combn(2) %>% 
  t %>% 
  data.frame %>% 
  rename(
    week1 = X1, 
    week2 = X2
  )

grds <- allind %>% 
  select(trt, species, var) %>% 
  unique() %>% 
  crossing(wks) %>% 
  mutate(
    res = purrr::pmap(list(trt, species, var, week1, week2), function(trt, species, var, week1, week2){

      specsel <- species
      varsel <- var
      trtsel <- trt
      
      tomod <- allind %>%
        filter(species %in% specsel) %>%
        filter(var %in% varsel) %>%
        filter(trt %in% trtsel)

      val1 <- tomod %>%
        filter(week == week1) %>%
        pull(val)

      val2 <- tomod %>%
        filter(week == week2) %>%
        pull(val)

      est <- t.test(val1, val2)

      perc <-  as.numeric(100 * diff(est$estimate) / est$estimate[2])
      pout <- p_ast(est$p.value)
      out <- data.frame(perc, pout)

      return(out)

    })
  ) %>% 
  unnest('res')

Size

wgtvars <- c('Final size (cm2)')
toplo1 <- allind %>% 
  filter(var %in% wgtvars) %>% 
  group_by(week, trt, species, var) %>% 
  summarise(
    aveval = mean(val, na.rm = T), 
    lowval = t.test(val)$conf.int[1],
    uppval = t.test(val)$conf.int[2]
  )

p1 <- ggplot(toplo1, aes(x = week, y = aveval)) + 
  geom_point(size = 2) + 
  geom_errorbar(aes(ymax = uppval, ymin = lowval), width = 1) + 
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 2, 4, 6)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Average values and 95% confidence intervals across time periods', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Average (+/- 95% CI)', 
    x = 'Week'
  )

toplo2 <- grds %>% 
  filter(var %in% wgtvars)

p2 <- ggplot(toplo2, aes(x = week1, y = week2, fill = perc)) + 
  geom_tile(color = 'black') + 
  geom_text(aes(label = pout)) + 
  scale_x_continuous('Week', expand = c(0, 0), breaks = c(0, 2, 4)) +
  scale_y_continuous('Week', expand = c(0, 0), breaks = c(2, 4, 6)) +
  # scale_fill_distiller('% diff.', palette = 'Purples', direction = 1) + 
  scale_fill_gradient2(low = 'lightgreen', mid = 'white', high = 'lightblue') +
  facet_grid(trt ~ species + var) +
  theme(
    panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank()
  ) + 
  labs(
    title = 'Differences in means between experiment time points', 
    subtitle = 'Significant differences are indicated by text'
  )

toplo3 <- allind %>% 
  filter(var %in% wgtvars) 

p3 <- ggplot(toplo3, aes(x = val)) + 
  stat_ecdf(aes(colour = trt), size = 1) +
  facet_grid(week ~ species + var, scales = 'free') +
  scale_colour_viridis_d() +
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank(), 
    legend.title = element_blank()
  ) + 
  labs(
    title = 'Cumulative distribution plots', 
    subtitle = 'Panels separated by species, week, and response measure',
    y = 'Cumulative frequency', 
    x = 'Weight (g)'
  )

p1 + p2 + p3 + plot_layout(ncol = 2)

Size across individuals

Observed

# all individuals
datprp <- read_csv(here::here('data/raw/finaldata-with-newmissingdata.csv')) %>% 
  clean_names() %>%
  select(week, trt = treatment, jar, id = individual_id, species, final_size = surface_area_cm2_final, initial_size = surface_area_cm2_initial)

toplo1 <- datprp %>% 
  gather('var', 'val', initial_size, final_size) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  mutate(
    week = case_when(
      period == 'initial' ~ 0, 
      T ~ week
    )
  ) %>% 
  unite('jar_id', jar, id)

p1 <- ggplot(toplo1, aes(x = week, y = val, group = jar_id)) + 
  geom_line() + 
  geom_point(size = 2) + 
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 2, 4, 6)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking size by individuals from time zero', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Size (cm2)', 
    x = 'Week'
  )

toplo2 <- datprp %>% 
  mutate(
    delt_size = final_size - initial_size, 
    init_size = 0,
  ) %>% 
  filter(week != 0) %>% 
  select(week, trt, jar, id, species, delt_size, init_size) %>% 
  gather('var', 'val', delt_size, init_size) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  filter(!week == 0) %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0, 
      T ~ week
    )
  ) %>% 
  unite('jar_id', jar, id)

p2 <- ggplot(toplo2, aes(x = week, y = val, group = jar_id)) + 
  geom_line() + 
  geom_point(size = 2) + 
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 2, 4, 6)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking size by individuals from time zero', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Size change from time zero (g)', 
    x = 'Week'
  ) + 
  geom_hline(yintercept = 0, linetype = 'dotted', colour = 'black')

toplo3 <- datprp %>% 
  mutate(
    scld_size = (final_size - initial_size) / week + initial_size,
    init_size = initial_size,
  ) %>% 
  filter(week != 0) %>% 
  select(week, trt, jar, id, species, scld_size, init_size) %>% 
  gather('var', 'val', scld_size, init_size) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  filter(!week == 0) %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0L, 
      T ~ 1L
    )
  ) %>% 
  unite('jar_id', jar, id)

p3 <- ggplot(toplo3, aes(x = week, y = val, group = jar_id)) + 
  geom_line(size = 0.5) +
  geom_point(size = 2) + 
  # geom_point(aes(y = aveval), colour = 'black', size = 3) +
  # geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  # geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled size by individuals from time zero, all weeks', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Scaled size (cm2/week)', 
    x = 'Week'
  ) 

toplo4 <- toplo3 %>% 
  group_by(week, trt, species, period, var) %>% 
  summarize(
    aveval = mean(val, na.rm = T), 
    uprval = t.test(val)$conf.int[2],
    lwrval = t.test(val)$conf.int[1]
  )

p4 <- ggplot(toplo4, aes(x = week)) + 
  geom_point(aes(y = aveval), colour = 'black', size = 3) +
  geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled size by individuals from time zero, all weeks',  
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Average scaled size (cm2/week) +/- 95% CI', 
    x = 'Week'
  ) 

toplo5 <- datprp %>% 
  mutate(
    scld_size = (final_size - initial_size) / week + initial_size,
    init_size = initial_size,
  ) %>% 
  filter(week == 2) %>% 
  select(week, trt, jar, id, species, scld_size, init_size) %>% 
  gather('var', 'val', scld_size, init_size) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0L, 
      T ~ 1L
    )
  ) %>% 
  unite('jar_id', jar, id)

p5 <- ggplot(toplo5, aes(x = week, y = val, group = jar_id)) + 
  geom_line(size = 0.5) +
  geom_point(size = 2) + 
  # geom_point(aes(y = aveval), colour = 'black', size = 3) +
  # geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  # geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled size by individuals from time zero, week two only', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Scaled size (cm2/week)', 
    x = 'Week'
  ) 

toplo6 <- toplo5 %>% 
  group_by(week, trt, species, period, var) %>% 
  summarize(
    aveval = mean(val, na.rm = T), 
    uprval = t.test(val)$conf.int[2],
    lwrval = t.test(val)$conf.int[1]
  )

p6 <- ggplot(toplo6, aes(x = week)) + 
  geom_point(aes(y = aveval), colour = 'black', size = 3) +
  geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled size by individuals from time zero, week two only',  
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Average Scaled size (cm2/week) +/- 95% CI', 
    x = 'Week'
  ) 

toplo7 <- datprp %>% 
  mutate(
    scld_size = (final_size - initial_size) / week + initial_size,
    init_size = initial_size,
  ) %>% 
  filter(week == 4) %>% 
  select(week, trt, jar, id, species, scld_size, init_size) %>% 
  gather('var', 'val', scld_size, init_size) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0L, 
      T ~ 1L
    )
  ) %>% 
  unite('jar_id', jar, id)

p7 <- ggplot(toplo7, aes(x = week, y = val, group = jar_id)) + 
  geom_line(size = 0.5) +
  geom_point(size = 2) + 
  # geom_point(aes(y = aveval), colour = 'black', size = 3) +
  # geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  # geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled size by individuals from time zero, week four only', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Scaled size (cm2/week)', 
    x = 'Week'
  ) 

toplo8 <- toplo7 %>% 
  group_by(week, trt, species, period, var) %>% 
  summarize(
    aveval = mean(val, na.rm = T), 
    uprval = t.test(val)$conf.int[2],
    lwrval = t.test(val)$conf.int[1]
  )

p8 <- ggplot(toplo8, aes(x = week)) + 
  geom_point(aes(y = aveval), colour = 'black', size = 3) +
  geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled size by individuals from time zero, week four only',  
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Average Scaled size (cm2/week) +/- 95% CI', 
    x = 'Week'
  ) 

toplo9 <- datprp %>% 
  mutate(
    scld_size = (final_size - initial_size) / week + initial_size,
    init_size = initial_size,
  ) %>% 
  filter(week == 6) %>% 
  select(week, trt, jar, id, species, scld_size, init_size) %>% 
  gather('var', 'val', scld_size, init_size) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0L, 
      T ~ 1L
    )
  ) %>% 
  unite('jar_id', jar, id)

p9 <- ggplot(toplo9, aes(x = week, y = val, group = jar_id)) + 
  geom_line(size = 0.5) +
  geom_point(size = 2) + 
  # geom_point(aes(y = aveval), colour = 'black', size = 3) +
  # geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  # geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled size by individuals from time zero, week six only', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Scaled size (cm2/week)', 
    x = 'Week'
  ) 

toplo10 <- toplo9 %>% 
  group_by(week, trt, species, period, var) %>% 
  summarize(
    aveval = mean(val, na.rm = T), 
    uprval = t.test(val)$conf.int[2],
    lwrval = t.test(val)$conf.int[1]
  )

p10 <- ggplot(toplo4, aes(x = week)) + 
  geom_point(aes(y = aveval), colour = 'black', size = 3) +
  geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled size by individuals from time zero, week six only',  
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Average Scaled size (cm2/week) +/- 95% CI', 
    x = 'Week'
  ) 

toplo11 <- datprp %>% 
  mutate(
    delt_size = final_size - initial_size, 
    init_size = 0,
  ) %>% 
  filter(week != 0) %>% 
  select(week, trt, jar, id, species, delt_size, init_size) %>% 
  gather('var', 'val', delt_size, init_size) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  filter(!week == 0) %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0, 
      T ~ week
    )
  ) %>% 
  unite('jar_id', jar, id) %>% 
  filter(week != 0) %>% 
  mutate(
    val = case_when(
      sign(val) %in% c(0, 1) ~ 1, 
      sign(val) == -1 ~ -1
    )
  ) %>% 
  group_by(week, trt, species, period, var, val) %>% 
  summarise(cnt = sum(val)) %>% 
  na.omit() %>% 
  ungroup %>% 
  mutate(val = factor(val, levels = c(1, -1), labels = c('positive', 'negative')))

p11 <- ggplot(toplo11, aes(x = week, y = cnt, fill = factor(val))) + 
  geom_bar(stat = 'identity') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(2, 4, 6)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank(), 
    legend.title = element_blank()
  ) + 
  labs(
    title = 'Tracking size by individuals from time zero', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Count of individuals with positive or negative size change', 
    x = 'Week'
  )

p1 + p2 + p11 + plot_spacer() + p3 + p4 + p5 + p6 + p7 + p8 + p9 + p10 + plot_layout(ncol = 2)

Modelled, all changes

Olympia oyster

indbymods %>% 
  filter(species == 'Olympia') %>% 
  filter(var == 'Final size (cm2)') %>% 
  pull(plomod)

Pacific oyster

indbymods %>% 
  filter(species == 'Pacific') %>% 
  filter(var == 'Final size (cm2)') %>% 
  pull(plomod)

Modelled, positive size

Olympia oyster

dirindbymods %>% 
  filter(species == 'Olympia') %>% 
  filter(var == 'Final size (cm2)') %>% 
  filter(direction == 1) %>% 
  pull(plomod)

Pacific oyster

dirindbymods %>% 
  filter(species == 'Pacific') %>% 
  filter(var == 'Final size (cm2)') %>% 
  filter(direction == 1) %>% 
  pull(plomod)

Modelled, negative size

Olympia oyster

dirindbymods %>% 
  filter(species == 'Olympia') %>% 
  filter(var == 'Final size (cm2)') %>% 
  filter(direction == -1) %>% 
  pull(plomod)

Pacific oyster

dirindbymods %>% 
  filter(species == 'Pacific') %>% 
  filter(var == 'Final size (cm2)') %>% 
  filter(direction == -1) %>% 
  pull(plomod)

Weight

wgtvars <- c('Shell weight (g)', 'Tissue weight (g)')
toplo1 <- allind %>% 
  filter(var %in% wgtvars) %>% 
  group_by(week, trt, species, var) %>% 
  summarise(
    aveval = mean(val, na.rm = T), 
    lowval = t.test(val)$conf.int[1],
    uppval = t.test(val)$conf.int[2]
  )

p1 <- ggplot(toplo1, aes(x = week, y = aveval)) + 
  geom_point(size = 2) + 
  geom_errorbar(aes(ymax = uppval, ymin = lowval), width = 1) + 
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 2, 4, 6)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Average values and 95% confidence intervals across time periods', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Average (+/- 95% CI)', 
    x = 'Week'
  )

toplo2 <- grds %>% 
  filter(var %in% wgtvars)

p2 <- ggplot(toplo2, aes(x = week1, y = week2, fill = perc)) + 
  geom_tile(color = 'black') + 
  geom_text(aes(label = pout)) + 
  scale_x_continuous('Week', expand = c(0, 0), breaks = c(0, 2, 4)) +
  scale_y_continuous('Week', expand = c(0, 0), breaks = c(2, 4, 6)) +
  # scale_fill_distiller('% diff.', palette = 'Purples', direction = 1) + 
  scale_fill_gradient2(low = 'lightgreen', mid = 'white', high = 'lightblue') +
  facet_grid(trt ~ species + var) +
  theme(
    panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank()
  ) + 
  labs(
    title = 'Differences in means between experiment time points', 
    subtitle = 'Significant differences are indicated by text'
  )

toplo3 <- allind %>% 
  filter(var %in% wgtvars) 

p3 <- ggplot(toplo3, aes(x = val)) + 
  stat_ecdf(aes(colour = trt), size = 1) +
  facet_grid(week ~ species + var, scales = 'free') +
  scale_colour_viridis_d() +
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank(), 
    legend.title = element_blank()
  ) + 
  labs(
    title = 'Cumulative distribution plots', 
    subtitle = 'Panels separated by species, week, and response measure',
    y = 'Cumulative frequency', 
    x = 'Weight (g)'
  )

p1 + p2 + p3 + plot_layout(ncol = 2)

Weight across individuals

Observed

# calc initial final from sum of tissue and shell weight
datprp <- read.csv(here::here('data/raw/Weight_with dead but not multiple_Kelp_4_3.csv'), stringsAsFactors = F) %>% 
  clean_names() %>% 
  # filter(week != 0) %>% 
  filter(species != '') %>% 
  mutate(
    species = gsub('\\s*$', '', species),
    whole_organism_weight = case_when(
      is.na(whole_organism_weight) &  !is.na(shell_weight) & !is.na(tissue_weight) ~ shell_weight + tissue_weight, 
      T ~ whole_organism_weight
    )
  ) %>% 
  filter(!is.na(initial_wetweight) & !is.na(whole_organism_weight)) %>% 
  select(week, trt = treatment, jar = i_jar, id = individual_id, species, initial_weight = initial_wetweight, 
         final_weight = whole_organism_weight) 

toplo1 <- datprp %>% 
  gather('var', 'val', initial_weight, final_weight) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  mutate(
    week = case_when(
      period == 'initial' ~ 0L, 
      T ~ week
    )
  ) %>% 
  unite('jar_id', jar, id)

p1 <- ggplot(toplo1, aes(x = week, y = val, group = jar_id)) + 
  geom_line() + 
  geom_point(size = 2) + 
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 2, 4, 6)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking weight by individuals from time zero', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Weight (g)', 
    x = 'Week'
  )

toplo2 <- datprp %>% 
  mutate(
    delt_weight = final_weight - initial_weight, 
    init_weight = 0,
  ) %>% 
  filter(week != 0) %>% 
  select(week, trt, jar, id, species, delt_weight, init_weight) %>% 
  gather('var', 'val', delt_weight, init_weight) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  filter(!week == 0) %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0L, 
      T ~ week
    )
  ) %>% 
  unite('jar_id', jar, id)

p2 <- ggplot(toplo2, aes(x = week, y = val, group = jar_id)) + 
  geom_line() + 
  geom_point(size = 2) + 
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 2, 4, 6)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking weight by individuals from time zero', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Weight change from time zero (g)', 
    x = 'Week'
  ) + 
  geom_hline(yintercept = 0, linetype = 'dotted', colour = 'black')

toplo3 <- datprp %>% 
  mutate(
    scld_weight = (final_weight - initial_weight) / week + initial_weight,
    init_weight = initial_weight,
  ) %>% 
  filter(week != 0) %>% 
  select(week, trt, jar, id, species, scld_weight, init_weight) %>% 
  gather('var', 'val', scld_weight, init_weight) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  filter(!week == 0) %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0L, 
      T ~ 1L
    )
  ) %>% 
  unite('jar_id', jar, id)

p3 <- ggplot(toplo3, aes(x = week, y = val, group = jar_id)) + 
  geom_line(size = 0.5) +
  geom_point(size = 2) + 
  # geom_point(aes(y = aveval), colour = 'black', size = 3) +
  # geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  # geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled weight by individuals from time zero, all weeks', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Scaled weight (g/week)', 
    x = 'Week'
  ) 

toplo4 <- toplo3 %>% 
  group_by(week, trt, species, period, var) %>% 
  summarize(
    aveval = mean(val, na.rm = T), 
    uprval = t.test(val)$conf.int[2],
    lwrval = t.test(val)$conf.int[1]
  )

p4 <- ggplot(toplo4, aes(x = week)) + 
  geom_point(aes(y = aveval), colour = 'black', size = 3) +
  geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled weight by individuals from time zero, all weeks',  
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Average scaled weight (g/week) +/- 95% CI', 
    x = 'Week'
  ) 

toplo5 <- datprp %>% 
  mutate(
    scld_weight = (final_weight - initial_weight) / week + initial_weight,
    init_weight = initial_weight,
  ) %>% 
  filter(week == 2) %>% 
  select(week, trt, jar, id, species, scld_weight, init_weight) %>% 
  gather('var', 'val', scld_weight, init_weight) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0L, 
      T ~ 1L
    )
  ) %>% 
  unite('jar_id', jar, id)

p5 <- ggplot(toplo5, aes(x = week, y = val, group = jar_id)) + 
  geom_line(size = 0.5) +
  geom_point(size = 2) + 
  # geom_point(aes(y = aveval), colour = 'black', size = 3) +
  # geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  # geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled weight by individuals from time zero, week two only', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Scaled weight (g/week)', 
    x = 'Week'
  ) 

toplo6 <- toplo5 %>% 
  group_by(week, trt, species, period, var) %>% 
  summarize(
    aveval = mean(val, na.rm = T), 
    uprval = t.test(val)$conf.int[2],
    lwrval = t.test(val)$conf.int[1]
  )

p6 <- ggplot(toplo6, aes(x = week)) + 
  geom_point(aes(y = aveval), colour = 'black', size = 3) +
  geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled weight by individuals from time zero, week two only',  
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Average Scaled weight (g/week) +/- 95% CI', 
    x = 'Week'
  ) 

toplo7 <- datprp %>% 
  mutate(
    scld_weight = (final_weight - initial_weight) / week + initial_weight,
    init_weight = initial_weight,
  ) %>% 
  filter(week == 4) %>% 
  select(week, trt, jar, id, species, scld_weight, init_weight) %>% 
  gather('var', 'val', scld_weight, init_weight) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0L, 
      T ~ 1L
    )
  ) %>% 
  unite('jar_id', jar, id)

p7 <- ggplot(toplo7, aes(x = week, y = val, group = jar_id)) + 
  geom_line(size = 0.5) +
  geom_point(size = 2) + 
  # geom_point(aes(y = aveval), colour = 'black', size = 3) +
  # geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  # geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled weight by individuals from time zero, week four only', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Scaled weight (g/week)', 
    x = 'Week'
  ) 

toplo8 <- toplo7 %>% 
  group_by(week, trt, species, period, var) %>% 
  summarize(
    aveval = mean(val, na.rm = T), 
    uprval = t.test(val)$conf.int[2],
    lwrval = t.test(val)$conf.int[1]
  )

p8 <- ggplot(toplo8, aes(x = week)) + 
  geom_point(aes(y = aveval), colour = 'black', size = 3) +
  geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled weight by individuals from time zero, week four only',  
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Average Scaled weight (g/week) +/- 95% CI', 
    x = 'Week'
  ) 

toplo9 <- datprp %>% 
  mutate(
    scld_weight = (final_weight - initial_weight) / week + initial_weight,
    init_weight = initial_weight,
  ) %>% 
  filter(week == 6) %>% 
  select(week, trt, jar, id, species, scld_weight, init_weight) %>% 
  gather('var', 'val', scld_weight, init_weight) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0L, 
      T ~ 1L
    )
  ) %>% 
  unite('jar_id', jar, id)

p9 <- ggplot(toplo9, aes(x = week, y = val, group = jar_id)) + 
  geom_line(size = 0.5) +
  geom_point(size = 2) + 
  # geom_point(aes(y = aveval), colour = 'black', size = 3) +
  # geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  # geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled weight by individuals from time zero, week six only', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Scaled weight (g/week)', 
    x = 'Week'
  ) 

toplo10 <- toplo9 %>% 
  group_by(week, trt, species, period, var) %>% 
  summarize(
    aveval = mean(val, na.rm = T), 
    uprval = t.test(val)$conf.int[2],
    lwrval = t.test(val)$conf.int[1]
  )

p10 <- ggplot(toplo4, aes(x = week)) + 
  geom_point(aes(y = aveval), colour = 'black', size = 3) +
  geom_errorbar(aes(ymin = lwrval, ymax = uprval), colour = 'black', width = 0.1) +
  geom_line(aes(y = aveval), colour = 'black') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(0, 1)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank()
  ) + 
  labs(
    title = 'Tracking scaled weight by individuals from time zero, week six only',  
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Average Scaled weight (g/week) +/- 95% CI', 
    x = 'Week'
  ) 

toplo11 <- datprp %>% 
  mutate(
    delt_weight = final_weight - initial_weight, 
    init_weight = 0,
  ) %>% 
  filter(week != 0) %>% 
  select(week, trt, jar, id, species, delt_weight, init_weight) %>% 
  gather('var', 'val', delt_weight, init_weight) %>% 
  separate(var, c('period', 'var')) %>% 
  group_by(week, trt, jar, id, species, period, var) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  filter(!week == 0) %>% 
  mutate(
    week = case_when(
      period == 'init' ~ 0L, 
      T ~ week
    )
  ) %>% 
  unite('jar_id', jar, id) %>% 
  filter(week != 0) %>% 
  mutate(
    val = case_when(
      sign(val) %in% c(0, 1) ~ 1, 
      sign(val) == -1 ~ -1
    )
  ) %>% 
  group_by(week, trt, species, period, var, val) %>% 
  summarise(cnt = sum(val)) %>% 
  na.omit() %>% 
  ungroup %>% 
  mutate(val = factor(val, levels = c(1, -1), labels = c('positive', 'negative')))

p11 <- ggplot(toplo11, aes(x = week, y = cnt, fill = factor(val))) + 
  geom_bar(stat = 'identity') +
  facet_grid(trt ~ species + var) +
  scale_x_continuous(breaks = c(2, 4, 6)) + 
  theme_bw() + 
  theme(
    # panel.background = element_blank(), 
    strip.background = element_blank(),
    axis.ticks = element_blank(), 
    panel.grid.minor = element_blank(), 
    legend.title = element_blank()
  ) + 
  labs(
    title = 'Tracking weight by individuals from time zero', 
    subtitle = 'Panels separated by species, treatment, and response measure',
    y = 'Count of individuals with positive or negative weight change', 
    x = 'Week'
  )

p1 + p2 + p11 + plot_spacer() + p3 + p4 + p5 + p6 + p7 + p8 + p9 + p10 + plot_layout(ncol = 2)

Modelled, all changes

Olympia oyster

indbymods %>% 
  filter(species == 'Olympia') %>% 
  filter(var == 'Whole weight (g)') %>% 
  pull(plomod)

Pacific oyster

indbymods %>% 
  filter(species == 'Pacific') %>% 
  filter(var == 'Whole weight (g)') %>% 
  pull(plomod)

Modelled, positive weight

Olympia oyster

dirindbymods %>% 
  filter(species == 'Olympia') %>% 
  filter(var == 'Whole weight (g)') %>% 
  filter(direction == 1) %>% 
  pull(plomod)

Pacific oyster

dirindbymods %>% 
  filter(species == 'Pacific') %>% 
  filter(var == 'Whole weight (g)') %>% 
  filter(direction == 1) %>% 
  pull(plomod)

Modelled, negative weight

Olympia oyster

dirindbymods %>% 
  filter(species == 'Olympia') %>% 
  filter(var == 'Whole weight (g)') %>% 
  filter(direction == -1) %>% 
  pull(plomod)

Pacific oyster

dirindbymods %>% 
  filter(species == 'Pacific') %>% 
  filter(var == 'Whole weight (g)') %>% 
  filter(direction == -1) %>% 
  pull(plomod)

Size vs weight vs. length

toplo1 <- allind %>% 
  group_by_at(vars(-val)) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  spread(var, val)

p1 <- ggplot(toplo1, aes(x = `Final size (cm2)`, y = `Shell weight (g)`, colour = trt)) +
  geom_point(alpha = 0.7) +
  # geom_smooth(se = F, method = 'lm') +
  scale_colour_viridis_d() +
  facet_grid(week ~ species) +
  theme_bw() +
  theme(
    # panel.background = element_blank(),
    strip.background = element_blank(),
    axis.ticks = element_blank(),
    panel.grid.minor = element_blank(),
    legend.title = element_blank()
  ) +
  labs(
    title = 'Tracking shell weight vs size',
    subtitle = 'Panels separated by species, week, and treatment',
    y = 'Shell weight (g)',
    x = 'Final size (cm2)'
  )

p2 <- ggplot(toplo1, aes(x = `Final size (cm2)`, y = `Tissue weight (g)`, colour = trt)) +
  geom_point(alpha = 0.7) +
  # geom_smooth(se = F, method = 'lm') +
  scale_colour_viridis_d() +
  facet_grid(week ~ species) +
  theme_bw() +
  theme(
    # panel.background = element_blank(),
    strip.background = element_blank(),
    axis.ticks = element_blank(),
    panel.grid.minor = element_blank(),
    legend.title = element_blank()
  ) +
  labs(
    title = 'Tracking tissue weight by size',
    subtitle = 'Panels separated by species, week, and treatment',
    y = 'Tissue weight (g)',
    x = 'Final size (cm2)'
  )

p3 <- ggplot(toplo1, aes(x = `Shell weight (g)`, y = `Tissue weight (g)`, colour = trt)) +
  geom_point(alpha = 0.7) +
  # geom_smooth(se = F, method = 'lm') +
  scale_colour_viridis_d() +
  facet_grid(week ~ species) +
  theme_bw() +
  theme(
    # panel.background = element_blank(),
    strip.background = element_blank(),
    axis.ticks = element_blank(),
    panel.grid.minor = element_blank(),
    legend.title = element_blank()
  ) +
  labs(
    title = 'Tracking tissue weight by shell weight',
    subtitle = 'Panels separated by species, week, and treatment',
    y = 'Tissue weight (g)',
    x = 'Shell weight (g)'
  )

lendat <- alllen <- read_excel(here::here('data/raw/Length_kelp_4_3_2020.xlsx')) %>%
  clean_names() %>%
  select(week, trt = treatment, jar, id = individual_id, species, final_length = length_cm_final, initial_length = length_cm_initial, final_width = width_cm_final, initial_width = width_cm_initial) %>%
  mutate(
    final_length = case_when(
      week == 0 & is.na(final_length) ~ initial_length,
      T ~ final_length
    ),
    final_width = case_when(
      week == 0 & is.na(final_width) ~ initial_width,
      T ~ final_width
    ),
    delt_length = final_length - initial_length,
    delt_width = final_width - initial_width,
    rate_length = case_when(
      week == 2 ~ delt_length / 14,
      week == 4 ~ delt_length / 28
    ),
    rate_width = case_when(
      week == 2 ~ delt_width / 14,
      week == 4 ~ delt_width / 28
    )
  ) %>%
  gather('var', 'val', -week, -trt, -jar, -id, -species) %>%
  filter(week != 6) %>%
  filter(!(week %in% 0 & var %in% c('delt_length', 'delt_width', 'rate_width', 'rate_length'))) %>%
  filter(!var %in% c('initial_length', 'initial_width')) %>% 
  group_by_at(vars(-val)) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  spread(var, val) %>% 
  mutate(jar = as.character(jar))

toplo2 <- allind %>% 
  group_by_at(vars(-val)) %>% 
  summarise(val = mean(val, na.rm = T)) %>% 
  ungroup %>% 
  spread(var, val) %>% 
  mutate(
    jar = as.character(jar), 
    trt = as.character(trt)
  ) %>% 
  inner_join(lendat, ., by = c('week', 'trt', 'jar', 'id', 'species')) %>% 
  mutate(
    trt = factor(trt, levels = trts$shrtlab, labels = trts$shrtlab)
  )

p4 <- ggplot(toplo2, aes(x = `Final size (cm2)`, y = final_length, colour = trt)) +
  geom_point(alpha = 0.7) +
  # geom_smooth(se = F, method = 'lm') +
  scale_colour_viridis_d() +
  facet_grid(week ~ species) +
  theme_bw() +
  theme(
    # panel.background = element_blank(),
    strip.background = element_blank(),
    axis.ticks = element_blank(),
    panel.grid.minor = element_blank(),
    legend.title = element_blank()
  ) +
  labs(
    title = 'Tracking length by size',
    subtitle = 'Panels separated by species, week, and treatment',
    y = 'Final length (cm)',
    x = 'Final size (cm2)'
  )

p5 <- ggplot(toplo2, aes(x = `Final size (cm2)`, y = final_width, colour = trt)) +
  geom_point(alpha = 0.7) +
  # geom_smooth(se = F, method = 'lm') +
  scale_colour_viridis_d() +
  facet_grid(week ~ species) +
  theme_bw() +
  theme(
    # panel.background = element_blank(),
    strip.background = element_blank(),
    axis.ticks = element_blank(),
    panel.grid.minor = element_blank(),
    legend.title = element_blank()
  ) +
  labs(
    title = 'Tracking width by size',
    subtitle = 'Panels separated by species, week, and treatment',
    y = 'Final width (cm)',
    x = 'Final size (cm2)'
  )

p6 <- ggplot(toplo2, aes(x = final_length, y = final_width, colour = trt)) +
  geom_point(alpha = 0.7) +
  # geom_smooth(se = F, method = 'lm') +
  scale_colour_viridis_d() +
  facet_grid(week ~ species) +
  theme_bw() +
  theme(
    # panel.background = element_blank(),
    strip.background = element_blank(),
    axis.ticks = element_blank(),
    panel.grid.minor = element_blank(),
    legend.title = element_blank()
  ) +
  labs(
    title = 'Tracking width by length',
    subtitle = 'Panels separated by species, week, and treatment',
    y = 'Final width (cm)',
    x = 'Final length (cm)'
  )

p1 + p2 + p3 + p4 + p5 + p6 + plot_layout(ncol = 2)

Survival

srvdat <- read.csv(here::here('data/raw/Weight_with dead but not multiple_Kelp_4_3.csv'), stringsAsFactors = F) %>% 
  clean_names() %>% 
  filter(species != '') %>% 
  mutate(species = gsub('\\s*$', '', species)) %>% 
  rename(
    id = individual_id, 
    jar = i_jar,
    trt = treatment
  ) %>% 
  select(week, trt, jar, id, species, dead) %>% 
  mutate(
    dead = case_when(
      dead %in% c('x', 'X') ~ 1, 
      T ~ 0
    ), 
    trt = factor(trt, levels = trts$shrtlab)
  ) %>% 
  group_by(species) %>%
  nest %>% 
  mutate(
    srvmod = purrr::map(data, function(x){
      
      fit <- survfit(Surv(week, dead) ~ trt, data = x)
      return(fit)
      
    }), 
    srvdif = purrr::map(data, function(x){
    
      fit <- survdiff(Surv(week, dead) ~ trt, data = x)
      return(fit)
      
    })
  )

srvplos <- srvdat %>% 
  mutate(
    plodat = purrr::map(srvmod, fortify)
  ) %>% 
  select(-data, -srvmod, -srvdif) %>% 
  unnest(plodat)

cols <- c('khaki4', 'khaki3', 'khaki1', 'seagreen4', 'seagreen2')
names(cols) <- trts$shrtlab

p <- ggplot(srvplos, aes(x = time, y = surv)) + 
  geom_line(aes(colour = strata), size = 1.5) + 
  geom_ribbon(aes(ymin = lower, ymax = upper), colour = NA, fill = 'lightgrey', alpha = 0.25) +
  facet_grid(species ~ strata) +
  scale_colour_manual('Treatment', values = cols, guide = F) +
  theme_minimal() +
  theme(
    panel.grid.major.x = element_blank(), 
    panel.grid.minor.x = element_blank(), 
    legend.position = 'top'
  ) + 
  labs(
    x = 'Week', 
    y = '% survival', 
    title = 'Survival estimates by treatment and species',
    subtitle = 'Shaded region is 95% CI'
  )
p

Olympia survival difference test

srvdat %>% 
  filter(species %in% 'Olympia') %>% 
  pull(srvdif) %>% 
  .[[1]]
## Call:
## survdiff(formula = Surv(week, dead) ~ trt, data = x)
## 
## n=1624, 6 observations deleted due to missingness.
## 
##               N Observed Expected (O-E)^2/E (O-E)^2/V
## trt=7.7C    323       30     27.0    0.3250    0.4300
## trt=7.7A0.2 324       21     27.0    1.3499    1.7860
## trt=7.7A0.5 329       29     26.8    0.1754    0.2316
## trt=8.0C    320       26     26.6    0.0125    0.0165
## trt=8.0A0.2 328       29     27.5    0.0800    0.1064
## 
##  Chisq= 2.1  on 4 degrees of freedom, p= 0.7

Pacific survival difference test

srvdat %>% 
  filter(species %in% 'Pacific') %>% 
  pull(srvdif) %>% 
  .[[1]]
## Call:
## survdiff(formula = Surv(week, dead) ~ trt, data = x)
## 
## n=1540, 5 observations deleted due to missingness.
## 
##               N Observed Expected (O-E)^2/E (O-E)^2/V
## trt=7.7C    303       22     20.9    0.0595    0.0770
## trt=7.7A0.2 301       21     22.0    0.0455    0.0596
## trt=7.7A0.5 319       21     22.6    0.1136    0.1498
## trt=8.0C    325       19     22.9    0.6541    0.8647
## trt=8.0A0.2 292       27     21.6    1.3253    1.7286
## 
##  Chisq= 2.3  on 4 degrees of freedom, p= 0.7